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ABSTRACT 

There  are  three  distinct  processes  by  which  upward-propagating  gravity  waves  influence  the  large-scale 
dynamics  and  energetics  of  the  middle  atmosphere:  (i)  nonlocalized  transport  of  momentum  through  wave 
propagation  in  three  dimensions  that  remotely  redistributes  atmospheric  momentum  in  both  zonal  and 
meridional  directions  from  wave  generation  to  wave  dissipation  regions;  (ii)  localized  diffusive  transport  of 
momentum,  heat,  and  tracers  due  to  mixing  induced  by  wave  breaking;  and  (iii)  localized  transport  of  heat  by 
perturbing  wave  structures  due  to  dissipation  that  redistributes  the  thermal  energy  within  a  finite  domain. 
These  effects  become  most  significant  for  breaking  waves  when  momentum  drag,  eddy  diffusion,  and  wave 
heating —  the  “breaking  trinity” — are  all  imposed  on  the  background  state.  This  paper  develops  a  3D  pa¬ 
rameterization  scheme  that  self-consistently  includes  the  breaking  trinity  in  large-scale  numerical  models. 
The  3D  parameterization  scheme  is  developed  based  on  the  general  relationship  between  the  wave  action  flux 
and  the  subgrid-scale  momentum  and  heat  fluxes  developed  by  Zhu  in  1987  and  a  mapping  approximation 
between  the  wave  source  spectrum  and  momentum  deposition  distribution  developed  by  Alexander  and 
Dunkerton  in  1999.  For  a  set  of  given  input  wind  and  temperature  profiles  at  each  model  grid,  the  parame¬ 
terization  scheme  outputs  the  vertical  profiles  of  the  subgrid-scale  force  terms  together  with  the  eddy  diffusion 
coefficients  in  the  momentum  and  energy  equations  for  a  3D  background  flow. 


1.  Introduction 

Three-dimensional  (3D)  numerical  models  are  impor¬ 
tant  tools  used  in  the  quantitative  study  of  the  middle 
atmosphere  and  interpretation  of  satellite  measurements. 
Well-designed  numerical  models  in  the  middle  atmo¬ 
sphere  often  consist  of  several  types  of  comprehensive 
modules  representing  radiation,  dynamics,  photochem¬ 
istry,  and  transport.  Among  these  four  types  of  modules, 
the  dynamics  or  the  dynamical  core  of  a  numerical  model 
plays  a  critical  role  in  organizing  and  coupling  different 
physical  processes  in  a  consistent  manner.  This  is  mainly 
due  to  the  dynamics  being  controlled  by  fluid  mechanics 
that  is  a  continuum,  behaving  highly  nonlinearly,  and 
evolving  with  time  on  many  different  scales.  The  other 
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modules  play  complementary  roles.  For  example,  the 
physics  module  for  calculating  radiative  heating  and 
photolysis  rates  provides  the  parameterized  force  to  drive 
the  dynamics  and  photochemistry.  This  partially  affects 
the  dynamics  with  bounded  or  controllable  uncertainty 
because  the  time-independent  problem  of  radiative 
transfer  is  well  defined  and  solvable  for  a  given  input  of 
atmospheric  parameters,  such  as  temperature  and  spe¬ 
cies,  solar  irradiance  spectrum,  and  kinetics  database 
(e.g.,  Zhu  2004;  Mlynczak  and  Zhou  1998;  Zhu  et  al. 
2007).  From  the  perspective  of  satellite  measurements 
and  the  global  modeling  of  atmospheric  structure,  it  is 
often  the  large-scale  dynamics  and  transport  that  can  be 
observed  directly  with  instruments  and  simulated  explic¬ 
itly  with  models.  For  example,  tidal  winds  were  directly 
observed  from  the  Upper  Atmosphere  Research  Satellite 
High  Resolution  Doppler  Imager  (UARS/HRDI)  and 
were  used  to  estimate  the  momentum  deposition  based 
on  the  derived  velocity  correlation  terms  that  character¬ 
ize  the  effect  of  tides  on  the  zonal  mean  flow  (Lieberman 
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and  Hays  1994).  The  temperature  measured  from  the 
Sounding  of  the  Atmosphere  using  Broadband  Emission 
Radiometer  (SABER)  onboard  the  Thermosphere,  Ion¬ 
osphere,  Mesosphere,  Energetics  and  Dynamics  (TIMED) 
satellite  has  been  used  to  derive  planetary-scale  waves  for 
both  temperature  and  winds  and  the  associated  wave  force 
terms  such  as  the  Eliassen-Palm  flux  divergence  of  tides 
in  the  middle  atmosphere  (Zhu  et  al.  2005,  2008).  These 
large-scale  fields  derived  from  satellite  observations  can 
be  directly  assimilated  into  numerical  models  in  a  di¬ 
agnostic  analysis  to  gain  additional  physical  insights  into 
middle  atmosphere  dynamics  (e.g.,  Akmaev  1997;  Zhu 
et  al.  2005).  On  the  other  hand,  the  momentum  and  en¬ 
ergy  sources  produced  by  subgrid-scale  motions  have  to 
be  parameterized. 

For  middle  atmospheric  modeling  studies,  the  subgrid- 
scale  motions  can  be  most  effectively  represented  by 
gravity  waves  because  of  their  ability  to  transport  mo¬ 
mentum  and  energy  over  a  large  spatial  distance.  Gen¬ 
erated  mainly  in  the  troposphere  by  mechanisms  such  as 
flow  over  topography  (e.g.,  Nappo  2002;  Teixeira  et  al. 
2004),  convection  (e.g.,  Alexander  et  al.  1995),  shear  in¬ 
stability  (e.g.,  Lindzen  1974;  Scinocca  and  Ford  2000),  and 
geostrophic  adjustment  (e.g.,  Zhu  and  Holton  1987), 
gravity  waves  can  propagate  upward  into  the  mesosphere 
and  lower  thermosphere  (MLT)  (e.g.,  Hines  1960).  It  has 
been  well  established  that  zonal  forces  generated  by 
breaking  gravity  waves  are  crucial  to  maintain  large-scale 
dynamics  and  transport  in  the  middle  atmosphere,  espe¬ 
cially  in  the  MLT  region  (e.g.,  Lindzen  1981;  Holton  1983; 
Fritts  1984;  Holton  and  Zhu  1984;  McIntyre  2000;  Holton 
and  Alexander  2000).  It  is  also  known,  both  observa- 
tionally  and  theoretically,  that  wave  generation  and 
propagation  are  3D  in  space  (e.g.,  Fritts  1984,  1995;  Zhu 
1987;  Sato  1994;  Marks  and  Eckermann  1995;  Broutman 
et  al.  2001,  2002,  2004).  On  a  slowly  varying  time  scale  or 
averaged  over  a  wave  period,  a  propagating  gravity  wave 
will  not  have  any  effect  on  the  background  wind  or  tem¬ 
perature  unless  the  wave  is  subject  to  dissipation  in  wave 
action  (e.g.,  Lighthill  1978;  Zhu  1987;  Eckermann  1992). 
On  the  other  hand,  a  dissipative  wave  could  (i)  induce 
a  momentum  drag  on  the  background  flow,  (ii)  produce 
eddy  mixing  in  the  background  dynamical  and  tracer 
fields,  and  (iii)  generate  a  heating-cooling  dipole  on  the 
background  temperature  field.  The  mechanisms  that  con¬ 
tribute  to  the  dissipation  of  a  gravity  wave  disturbance  in 
the  middle  atmosphere  are  convective  and  dynamical 
instabilities,  radiative  damping,  eddy  diffusion,  and 
nonlinear  wave-wave  interactions.  In  the  lower  ther¬ 
mosphere,  molecular  diffusion  makes  an  additional 
contribution.  Among  these  dissipative  processes,  the 
convective  and  dynamical  instabilities  associated  with 
wave  breaking  are  believed  to  be  dominant  in  producing 


momentum  drag  and  eddy  diffusion  in  the  background 
flow  (Fritts  and  Rastogi  1985).  This  is  mainly  because  the 
characteristic  gravity  waves  in  the  middle  atmosphere  are 
of  high  frequency  and  long  vertical  wavelength  (e.g., 
Hirota  and  Niki  1985;  Zhu  et  al.  1997),  such  that  they 
can  propagate  over  a  large  distance  without  significant 
damping  before  wave  breaking.  In  addition  to  the  drag 
and  eddy  diffusion  that  have  been  extensively  studied  and 
have  been  parameterized  in  many  middle  atmosphere 
models,  it  has  been  recognized  recently  that  the  dy¬ 
namical  heating  associated  with  the  downward  heat 
flux  by  dissipative  gravity  waves  also  needs  to  be  in¬ 
cluded  in  middle  atmosphere  models  (Becker  2004; 
Akmaev  2007).  Based  on  the  rocket  measurements  of 
neutral  density  fluctuations,  Liibken  (1997)  found  that 
the  dynamical  heating  rate  near  the  mesopause  estimated 
from  the  turbulence  energy  dissipation  rates  amounts 
to  10-20  K  day-1,  which  is  comparable  to  the  magni¬ 
tudes  of  radiative  and  chemical  heating  rates  in  the 
same  region.  In  summary,  the  most  effective  dissipa¬ 
tion  of  gravity  waves — wave  breaking — produces  three 
important  effects  on  the  background  state:  wave  drag, 
eddy  diffusion,  and  wave  heating.  Here,  we  call  these 
three  dynamical  and  energetic  consequences  resulting 
from  gravity  wave  breaking  in  the  middle  atmosphere 
the  “breaking  trinity.” 

In  terms  of  momentum  flux  or  Reynolds  stresses,  the 
parameterization  of  subgrid-scale  motions  by  gravity 
waves  is  fundamentally  different  from  those  in  so-called 
turbulent  viscosity  models,  where  the  Reynolds  stresses 
are  essentially  parameterized  by  various  length  theories, 
including  Prandtl’s  mixing-length  theory  and  iGepsilon 
models  (e.g.,  Tennekes  and  Lumley  1972;  Pope  2000).  In 
the  latter  case,  the  momentum  is  carried  and  transferred 
by  moving  fluid  parcels  for  a  localized  exchange.  For 
mechanistic  general  circulation  models  (GCMs)  with  a 
high  resolution  that  can  explicitly  resolve  medium-scale 
gravity  waves,  the  localized  mixing-length  theory  can 
also  be  adopted  to  parameterize  the  effects  of  eddies 
with  much  smaller  subgrid  scales  (Becker  2009).  On  the 
other  hand,  gravity  waves  carry  and  transfer  momentum 
by  pressure  and  velocity  fluctuations  over  a  large  spatial 
distance.  Such  a  nonlocalized  transport  of  momentum 
also  provides  one  of  the  important  physical  bases  for 
changes  in  the  tropospheric  circulation  to  affect  the 
MLT  dynamics  and  climatology.  From  the  perspective 
of  wave-mean  flow  interactions  in  3D  flow  on  how  the 
momentum  and  energy  are  spatially  redistributed  within 
the  entire  domain,  the  breaking  trinity  can  be  un¬ 
derstood  as  three  effects  of  wave  propagation  and  dis¬ 
sipation  on  the  3D  mean  flow:  (i)  nonlocalized  transport 
of  momentum  through  wave  propagation  in  3D  that 
remotely  redistributes  atmospheric  momentum  in  both 
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zonal  and  meridional  directions  from  wave  generation 
to  wave  dissipation  regions;  (ii)  localized  diffusive 
transport  of  momentum,  heat,  and  tracers  due  to  3D 
mixing  induced  by  wave  breaking;  and  (iii)  localized 
transport  of  heat  by  perturbing  wave  structures  that 
redistributes  the  thermal  energy  within  a  finite  domain. 
To  quantitatively  and  self-consistently  examine  these 
effects  on  the  3D  background  state,  we  need  a  parame¬ 
terization  scheme  that  can  provide  (i)  wave  drag  to  the 
background  3D  flow  consisting  of  both  the  zonal  and 
meridional  components  in  the  momentum  equations, 
(ii)  eddy  diffusion  to  the  dynamical  and  tracer  fields  that 
represents  the  mixing  or  the  field  smoothness  by  param¬ 
eterized  subgrid-scale  motions,  and  (iii)  wave  heating  to 
the  temperature  field  induced  by  the  wave  breaking  in  the 
energy  equation,  all  resulting  from  the  dissipation  of  wave 
action  caused  by  the  interaction  between  3D  subgrid-scale 
waves  and  3D  background  flow.  This  paper  introduces 
a  new  parameterization  scheme  of  drag,  eddy  diffusion, 
and  wave  heating  imposed  on  a  3D  background  state  by 
the  breaking  of  gravity  waves.  The  output  of  the  new 
parameterization  scheme  will  self-consistently  include 
all  three  distinct  processes  associated  with  the  breaking 
trinity  by  which  the  upward-propagating  gravity  waves 
influence  the  3D  large-scale  dynamics  and  energetics  of 
the  middle  atmosphere. 

In  section  2,  we  introduce  a  spectral  parameterization 
scheme  of  gravity  wave  breaking  that  simultaneously 
produces  the  breaking  trinity  for  given  3D  wind  and 
temperature  profiles.  The  central  part  of  the  derivation  is 
built  on  previous  work  by  Holton  and  Zhu  (1984),  Zhu 
(1987),  and  Alexander  and  Dunkerton  (1999),  which  was 
based  on  the  general  relationship  between  the  wave  ac¬ 
tion  flux  and  the  wave  momentum,  and  also  on  a  simple 
mapping  approximation  between  the  wave  source  spec¬ 
trum  in  wave  parameters  and  momentum  deposition 
distribution  in  altitude.  Section  3  shows  the  typical  mag¬ 
nitudes  and  patterns  of  the  breaking  trinity  to  a  3D 
background  field  of  wind  and  temperature  derived  from 
an  output  field  of  a  high-altitude  version  of  the  Goddard 
Earth  Observing  System  atmospheric  model  (GEOS-5). 
Finally,  section  4  summarizes  the  paper. 


which  can  be  resolved  explicitly  by  the  model  grid,  and 
the  other  is  for  3D  wave  fields.  For  simplicity  and  also 
because  the  effects  of  the  breaking  trinity  due  to  the 
subgrid-scale  motions  will  be  parameterized  at  each  model 
grid  independently,  we  develop  our  parameterization 
scheme  by  assuming  the  hydrostatic  approximation  on 
an  /  plane.  Using  the  two-scale  analysis  method,  the 
equations  for  the  evolution  of  a  background  state  under 
a  localized  wave  forcing  representing  the  subgrid-scale 
motions  can  be  written  as  (Holton  1975;  Zhu  1987) 
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In  the  above  equations,  u,  E,  and  w  are  the  background 
zonal  x,  meridional  y,  and  vertical  z  velocities,  re¬ 
spectively;  O  is  the  background  geopotential  and  is  hy¬ 
drostatically  related  to  the  background  temperature  T ; 
/is  the  Coriolis  parameter;  p0  is  the  background  density 
of  the  atmosphere;  N 2  is  the  squared  buoyancy  frequency; 
and  D  is  the  time  derivative  following  the  horizontal 
component  of  the  background  wind: 


— —  d  _  d  _  d 

D - +  u  —  +  v  — . 

dt  dx  dy 


(5) 


The  squared  buoyancy  frequency  slowly  varies  with  al¬ 
titude  and  is  a  measure  of  the  static  stability  of  the 
background  atmosphere: 


2.  Parameterization  of  gravity  wave  breaking  in 
a  three-dimensional  background  flow 

a.  Wave-mean  flow  interaction  in 
a  three-dimensional  atmosphere 

Following  Holton  (1975)  and  Zhu  (1987),  the  gravity 
wave-mean  flow  interactions  in  a  3D  atmosphere  are 
described  by  two  sets  of  coupled  equations.  One  is  for 
the  large-scale  and  slowly  varying  background  flow, 


N2(z ) 


RfdT 

H\Jz 


(6) 


where  R  is  the  gas  constant,  H  is  the  scale  height  of  the 
background  air  density,  k  =  R/cp  ~  2/7,  and  cp  is  the 
specific  heat  at  constant  pressure.  The  left-hand  sides  of 
Eqs.  (l)-(3)  represent  the  dynamical  evolution  of  a  flow 
on  a  grid-resolved  scale.  The  dynamical  cores  of  differ¬ 
ent  GCMs  may  differ  significantly.  The  right-hand  terms 
in  Eqs.  (l)-(3)  represent  the  momentum  and  energy 
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sources  contributed  by  the  nonconservative  processes  of 
dissipation  and  transience  of  upward-propagating  subgrid- 
scale  gravity  waves,  (SM,  Sv )  and  ST  are  the  wave  drag  and 
heating  terms,  and  Kzz-m  and  Kzz-T  are  the  eddy  dif¬ 
fusion  coefficients  for  the  momentum  and  sensible  heat, 
respectively.  Note  that  Kzz-m  and  Kzz-T  are  related  by 
the  eddy  Prandtl  number  Pr  (=  Kzz-m/Kzz-T).  Since  the 
sensible  heat  flux  is  proportional  to  the  potential  tem¬ 
perature  gradient,  Kzz-T  can  also  be  used  for  parameter¬ 
izing  subgrid-scale  chemical  tracer  transport.  In  Eq.  (4), 
(Fm,  ¥v)  and  are  the  pseudomomentum  flux  and  sen¬ 
sible  heat  flux,  respectively.  In  this  paper,  all  these  terms 
and  coefficients  are  calculated  for  a  spectral  distribution  of 
a  set  of  upward-propagating  gravity  waves  specified  at 
the  wave  source  region  in  the  lower  atmosphere. 
Computationally,  a  subroutine  that  uses  the  vertical 
wind  and  temperature  profiles  as  input  at  each  model 
grid  returns  five  vertical  profiles  of  Su,  Sv ,  ST ,  Kzz-m , 
and  Kzz-T  in  Eqs.  (l)-(3).  In  addition,  the  geographic 
location  and  seasonal  timing  also  determine  the  sta¬ 
tistical  properties  of  the  source  spectrum  for  the  gravity 
waves  specified  in  the  lower  atmosphere. 

In  Eqs.  (l)-(3),  the  momentum  (p0w,  p0v )  and  total 
potential  energy  (p0cpT)  under  the  Boussinesq  approx¬ 
imation  are  all  linear  variables,  so  the  eddy  forcing  terms 
on  the  right-hand  sides  of  the  mean  momentum  and 
energy  equations  can  all  be  expressed  in  flux  form,  as 
shown  in  Eq.  (4)  (Holton  1975).  As  a  result,  the  mean 
momentum  and  potential  energies  can  be  spatially  redis¬ 
tributed  through  the  flux  divergences,  but  the  integrations 
over  an  enclosed  domain  will  be  solely  determined  by  the 
boundary  fluxes.  Mathematically,  this  is  also  because  the 
average  of  a  perturbed  linear  quantity  vanishes  (e.g., 
pQu'  =  0)  so  the  eddy  terms  do  not  contribute  to  the  mean 
quantities  in  the  volume  integrations.  Some  parameteri¬ 
zation  schemes  for  gravity  wave  breaking  have  been  de¬ 
veloped  based  on  an  energy  equation  that  also  includes 
the  kinetic  energy  (e.g.,  Becker  2004).  Since  the  kinetic 
energy  is  a  quadratic  quantity,  the  separation  of  the  at¬ 
mospheric  motion  between  the  mean  and  eddies  will  re¬ 
sult  in  two  nonzero  terms  in  its  kinetic  energy  (e.g., 
pQu2  =  p0u2  +  p0u'2).  In  this  case,  the  mean  energy 
equation  can  no  longer  be  written  in  a  purely  flux  form 
because  there  exists  an  energy  conversion  between  the 
mean  and  eddy  fields.  Furthermore,  when  the  kinetic 
energy  is  to  be  explicitly  included  in  the  energy  budget, 
a  more  sensible  approach  is  to  introduce  the  available 
potential  energy  that  is  also  quadratic  in  temperature 
(Lorenz  1955).  Under  such  a  circumstance,  only  the 
total  energy  that  includes  both  the  mean  plus  eddy  and 
kinetic  plus  available  potential  energies  is  still  conserved 
and  can  be  written  in  a  purely  flux  form  (e.g.,  Dutton 
1976;  Holton  1975). 


The  effect  of  the  background  state  on  a  wave  com¬ 
ponent  is  described  by  a  set  of  ray  tracing  equations  in  a 
3D  stratified  flow  for  its  wave  parameters  (e.g.,  Lighthill 
1978;  Zhu  1987;  Marks  and  Eckermann  1995;  Broutman 
et  al.  2004).  Specifically,  the  wave  amplitude  is  best  de¬ 
scribed  by  the  wave  action  conservation  law 

—  +  V  •  (c  A)  =  —a  A,  (7) 

v  8  '  w  ’  v  7 

where  A  =  E/co  is  the  wave  action  density,  with  E  and  go 
being  the  wave  energy  density  and  intrinsic  frequency, 
respectively;  cg  is  the  group  velocity  and  c gA  is  called  the 
wave  action  flux;  and  aw  is  the  bulk  damping  rate  co¬ 
efficient  of  the  wave  component,  which  characterizes 
various  damping  processes  including  radiative  and  eddy 
diffusive  damping,  molecular  diffusion,  and  wave  break¬ 
ing.  A  wave  component  is  also  called  a  wave  packet  if 
its  wave  parameters  vary  slowly  with  space  and  time. 
The  interaction  between  a  wave  packet  and  the  back¬ 
ground  mean  flow  is  often  expressed  by  a  relationship 
between  the  wave  action  flux  and  pseudomomentum  flux 
or  Eliassen-Palm  flux  (e.g.,  Dunkerton  1981;  Palmer 
1982;  Andrews  et  al.  1987;  Alexander  and  Dunkerton 
1999).  Given  such  a  relationship,  one  can  calculate  the 
subgrid-scale  drag  terms  in  Eq.  (4)  once  the  wave  action 
density  has  been  solved  from  the  wave  action  conserva¬ 
tion  law  (7).  Since  the  drag  terms  and  dynamical  heating 
term  in  Eqs.  (l)-(4)  have  been  expressed  in  flux  form, 
the  conservations  of  momentum  and  thermal  energy  of 
the  dynamical  core  are  automatically  preserved  for  the 
grid-resolved  flow  as  long  as  no  spurious  sources  are 
introduced  on  the  boundaries.  Most  models  assume  a 
hydrostatic  approximation  for  both  large-scale  back¬ 
ground  state  and  subgrid-scale  wave  fields.  Recently, 
Shaw  and  Shepherd  (2009)  proposed  a  framework  to 
couple  a  hydrostatic  large-scale  flow  to  a  nonhydrostatic 
subgrid-scale  flow  that  preserves  the  momentum  and  to¬ 
tal  energy  conservations  self-consistently. 

Assuming  the  background  state  and  wave  parameters 
are  slowly  varying  in  time  and  space  with  small  dissi¬ 
pation,  neglecting  the  effect  of  planetary  rotation  on 
medium-scale  gravity  waves,  and  using  the  two-scale 
analysis  method,  Zhu  (1987)  showed  that  the  wave  ac¬ 
tion  density  and  the  pseudomomentum  and  sensible 
heat  fluxes  in  the  first-order  approximation  for  an  in¬ 
dividual  gravity  wave  packet  in  Eqs.  (l)-(4)  were  related 
through 

F«  =  kcgA  ’  F  =  lcoA  ’  f,6.  =  0,  (8a,b,c) 

where  k  and  /  are  the  horizontal  wavenumbers  of  the 
wave  packet  in  the  x  and  y  directions,  respectively. 
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-  ( u ,  v,  0)  is  the  intrinsic  group  velocity  Doppler- 
shifted  by  the  horizontal  component  of  the  background 
wind.  Equation  (8)  together  with  Eqs.  (l)-(4)  describe 
the  effect  of  a  wave  packet  on  the  background  state.  Note 
that  since  we  have  applied  the  hydrostatic  and  the  /-plane 
approximations  to  the  original  set  of  equations  there  is 
no  velocity  shift  in  the  vertical  component  of  the  group 
velocity,  although  it  varies  with  altitude:  c  =  c  (z). 
The  wave  dispersion  relation  and  cgz  under  the  hydro¬ 
static  approximation  are  given  by  (Zhu  1987) 


cb2  =  f2  + 


k2n2 


A2  +  0.25 !H2  ’ 


\K2N2 

Cgz  ~  _  &»(A2  +  0.25///2)’ 


(9) 

(10) 


where  A  is  the  vertical  wavenumber  and  K  is  the  total 
horizontal  wavenumber  defined  by  K2  =  k2  +  l2.  The 
inertial  effect  of  the  Coriolis  parameter  /  in  Eq.  (9) 
becomes  important  only  if  one  focuses  on  large-scale 
gravity  waves  with  horizontal  wavelengths  of  several 
hundred  kilometers  and  examines  their  effects  on  lateral 
propagations  (e.g.,  Dunkerton  1984).  Since  our  param¬ 
eterization  scheme  is  to  be  implemented  on  each  model 
grid  independently,  the  inertial  effect  is  often  neglected 
in  the  current  paper  except  when  its  effect  becomes 
important  in  prescribing  wave  parameters.  Also  note 
that  the  sensible  heat  flux  vanishes  under  the  first-order 
approximation  for  a  conservative  gravity  wave  compo¬ 
nent,  as  shown  above  in  Eq.  (8c).  Equation  (8)  shows  that 
for  a  given  wave  action  flux  c gA  that  has  been  shifted  by 
the  background  wind  the  ratio  of  the  momentum  fluxes 
between  x  and  y  directions  is  related  by  the  direction  of 
the  wavenumber  vector  (A,  /)  and  is  independent  of  the 
direction  of  the  background  flow  vector  [u(z),  u(z)].  This 
means  that  the  partition  of  the  pseudomomentum  flux 
vector  between  x  and  y  directions  is  determined  by  the 
wave  structure.  The  same  dependent  relations  of  the 
momentum  flux  on  the  wavenumber  vector  and  wave 
action  density  were  also  derived  by  Warner  and  McIntyre 
(1996)  in  their  extension  of  Lindzen’s  (1981)  parameter¬ 
ization  scheme,  where  the  distribution  and  variation  of 
the  energy  spectrum  were  prescribed  in  terms  of  A.  As¬ 
suming  a  steady  vertical  wind  profile  that  is  slowly  varying 
in  the  horizontal  direction,  it  can  be  shown  that  the  hor¬ 
izontal  wavenumber  (A,  /)  and  the  ground-relative  wave 
frequency  co  do  not  vary  with  altitude  (e.g.,  Zhu  1987; 
Eckermann  1992).  Therefore,  the  partition  between  ¥u 
and  ¥v  of  a  wave  component  is  entirely  determined  by  its 
horizontal  wave  vector  specified  at  the  source  region. 
Furthermore,  since  the  horizontal  components  of  the 


phase  velocity  (cx,  cy )  =  ( colk ,  co/l)  do  not  vary  with  alti¬ 
tude  either,  most  parameterization  schemes  have  speci¬ 
fied  the  source  spectra  as  a  function  of  (A,  /)  and  phase 
speed  c  at  the  lower  boundary  (e.g.,  Holton  and  Zhu  1984; 
Alexander  and  Dunkerton  1999;  Becker  and  Schmitz 
2002).  Such  a  dependence  on  wave  parameters  in  pre¬ 
scribing  the  wave  spectrum  at  the  source  region  will  also 
be  adopted  in  the  current  parameterization  scheme. 

The  wave  action  conservation  law  (7)  suggests  that  the 
wave  action  flux  (c gA)  is  a  conserved  quantity  in  the 
absence  of  wave  dissipation  (aw  =  0).  As  a  result,  A  is 
invariant  and  its  value  along  the  rays  depicted  by  cg  is 
a  constant  regardless  of  how  other  individual  wave  pa¬ 
rameters  vary  as  the  wave  packet  interacts  with  the 
background  state.  Specifically,  assuming  the  steady- 
state  approximation  (d/dt  =  0),  the  wave  action  flux  c gA 
at  two  different  cross  sections  in  space  will  be  identical 
and  thus  can  be  evaluated  easily.  On  the  other  hand,  the 
pseudomomentum  fluxes  in  Eq.  (8)  associated  with  the 
drag  terms  in  Eqs.  (1)  and  (2)  contain  the  intrinsic  group 
velocity,  which  prevents  one  from  directly  applying  the 
wave  action  conservation  law.  However,  for  a  parame¬ 
terization  scheme  in  the  middle  atmosphere  that  assumes 
the  upward-propagating  gravity  waves  are  dominant,  the 
momentum  flux  convergence  is  exclusively  contributed 
from  the  vertical  component  of  the  pseudomomentum 
fluxes.  Under  such  a  circumstance,  we  have 

F  =  cos OJKc  A ),  F  =  sin OJKc  A ), 

u,z  Cu  gz  v,z  (A  gz 

F ^  z  =  —com A,  (lla,b,c) 

where  we  have  introduced  the  horizontal  direction  of 
wave  propagation  60  based  on  the  horizontal  wave- 
number  vector  by  (A,  /)  =  K(cos00,  sin 0O).  We  have  also 
shown  the  higher-order  approximation  of  the  sensible 
heat  flux  in  Eq.  (lie),  where  co  and  cot  are  the  real  and 
imaginary  parts  of  the  intrinsic  frequency,  respectively 
(Holton  and  Zhu  1984).  Equation  (11)  shows  explicitly 
that  for  an  upward-propagating  ( cgz  >  0)  gravity  wave 
packet  the  signs  of  the  momentum  flux  components 
are  determined  by  the  wavenumber  vector.  Also,  the 
dissipation  ( cot  >  0)  induces  a  downward  heat  flux 
(F^  z  <  0)  for  an  upward-propagating  gravity  wave  that 
hash  positive  co  by  the  notation  convention  (Holton  and 
Zhu  1984). 

Expressions  similar  to  Eq.  (lie)  have  also  been  derived 
previously  (Walterscheid  1981;  Talaat  et  al.  2001;  Becker 
2004).  Here,  for  3D  applications  that  include  both  the 
zonal  and  meridional  velocity  components  and  also  for 
using  horizontal  wavenumber  and  phase  speed  as  work¬ 
ing  wave  parameters,  the  intrinsic  frequency  is  related  to 
the  ground-relative  frequency  co  by  co  =  co  —  ku  —  Iv.  If 
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we  further  define  the  direction  of  the  background  flow 
(u,v)  by  U(cos8 ,  sin5),  then  co  can  be  simplified  as 

co  =  cd-KU  cos  6  =  K(c  —  U  cos0),  (12) 

where  U  =  Vu2  +  v2  is  the  magnitude  of  the  back¬ 
ground  flow  and  0  =  60  —  8  is  the  angle  between  the 
horizontal  wavenumber  vector  and  background  flow. 
Again,  the  ground-relative  phase  speed  c  =  co/K  (>0) 
does  not  vary  with  altitude.  The  vertical  sensible  heat 
flux  by  the  gravity  wave  component  vanishes  unless  it  is 
subject  to  damping  or  wave  breaking  that  can  be  char¬ 
acterized  by  the  imaginary  part  of  its  frequency  w*.  Note 
that  it  is  the  flux  divergence  that  will  contribute  to  the 
momentum  and  energy  budgets  in  Eqs.  (l)-(3).  Fur¬ 
thermore,  the  flux  form  of  Eq.  (7)  suggests  that  the  in¬ 
ternal  damping  and  boundary  values  of  wave  action 
density  are  closely  related.  Hence,  wave  momentum  and 
heat  fluxes  make  similar  types  of  contributions  to  the 
momentum  and  energy  budgets  in  such  a  way  that  they 
only  redistribute  spatially  within  and  vanish  outside  the 
dissipation  and  generation  regions.  Specifically,  when 
the  heat  flux  for  a  wave  packet  F ^  z  vanishes  at  both 
the  lower  boundary  where  waves  are  conserved  (coi  =  0) 
and  the  upper  boundary  where  all  waves  have  been 
completely  dissipated  (A  =  0),  nonzero  wave  heating 
should  occur  in  the  form  of  a  dipole  so  that  the  wave 
heating  within  a  finite  domain  will  be  exactly  balanced 
by  wave  cooling  in  an  adjacent  region.  The  wave  damp¬ 
ing  associated  with  the  breaking  processes  also  leads 
to  turbulent  eddy  diffusion  that  can  be  parameterized 
as  eddy  diffusion  coefficients  for  the  background  state. 
Our  parameterization  scheme  provides  the  breaking 
trinity,  which  includes  both  the  momentum  drag  and 
sensible  heating  terms  by  wave  breaking  together  with 
the  eddy  diffusion  coefficients  Kzz-m  and  Kzz-T used  in 
Eqs.  (l)-(3). 

b.  Extension  of  Lindzen’s  parameterization  to 
a  three-dimensional  atmosphere 

We  have  already  indicated  that  the  convective  and 
dynamical  instabilities  associated  with  wave  breaking 
are  believed  to  be  the  dominant  mechanism  of  wave 
dissipation  for  the  upward-propagating  gravity  waves  in 
the  middle  atmosphere.  The  wave  breaking  can  produce 
both  the  momentum  drag  and  eddy  diffusion  to  the 
background  flow  (Lindzen  1981).  On  the  other  hand,  for 
slowly  propagating  planetary  waves  where  radiative 
damping  is  the  main  dissipation  mechanism,  the  effect  of 
eddy  diffusion  by  wave  dissipation  on  the  background 
state  is  small  and  often  neglected  (e.g.,  Matsuno  1970; 
Holton  and  Lindzen  1972).  Since  cot  represents  the  wave 
dissipation  effect  by  either  radiative  damping  or  eddy 


diffusion  (Holton  and  Zhu  1984),  the  wave  dynamical 
heating  associated  with  Eq.  (lie)  remains  potentially 
significant  for  a  slowly  propagating  wave  with  radiative 
damping  as  its  main  dissipation  mechanism.  The  in¬ 
stabilities  of  an  upward-propagating  gravity  wave  are 
mainly  caused  by  the  exponential  growth  of  wave  ampli¬ 
tude  due  to  the  atmospheric  density  effect  or  the  increase 
of  A  due  to  the  critical  level  and  are  modified  by  the  static 
stability  of  the  background  state.  The  most  widely  used 
parameterization  schemes  have  been  those  based  on 
Lindzen’s  (1981)  theory  of  wave  breaking  and  saturation. 
The  term  “saturation”  refers  to  the  growth  cessation  of  an 
upward-propagating  unstable  wave  component  when  its 
Richardson  number  reaches  its  critical  value  of  V4  (shear 
instability)  or  approaches  zero  (convective  instability). 
Lindzen’s  parameterization  scheme  (1981)  was  origi¬ 
nally  proposed  for  deriving  Su  and  the  associated  Kzz-m 
for  a  zonal  mean  flow  (w,  0).  Although  Lindzen’s  (1981) 
parameterization  scheme  was  extended  by  Holton  and 
Zhu  (1984)  and  several  others,  most  schemes  only  cal¬ 
culate  four  terms  including  Su,  Kzz-m ,  and  Kzz-T  in 
Eqs.  (l)-(3)  for  a  two-dimensional  (2D)  flow  such  as  the 
one  containing  the  x  component  of  the  wind  profile.  The 
central  aspect  in  Lindzen’s  (1981)  parameterization  is 
that  (i)  the  specification  of  the  momentum  deposition  of 
the  breaking  wave  component  between  the  breaking 
level  Zb  and  the  critical  level  zc  is  based  on  the  saturation 
assumption.  This  original  parameterization  scheme  was 
extended  by  Holton  (1982)  and  Holton  and  Zhu  (1984) 
in  the  following  two  aspects:  (ii)  a  directional  isotropic 
source  spectrum  of  waves  in  phase  speed  is  specified  at 
the  lower  boundary  and  (iii)  each  wave  component  is 
influenced  by  Newtonian  cooling  and  the  eddy  diffu¬ 
sion  induced  by  the  breaking  of  other  wave  compo¬ 
nents  with  lower  Zb-  Aspect  (iii)  allows  the  nonlinear 
interaction  among  different  wave  components  to  be 
partially  included  in  the  parameterization.  These  exten¬ 
sions  also  cause  the  wave  momentum  to  be  deposited  in 
a  more  extended  region  of  the  atmosphere,  especially 
below  the  dominant  wave  breaking  level  It  is  noted 
that  aspects  (i)  and  (iii)  critically  depend  on  aspect  (ii), 
which  varies  with  season  and  geography.  In  fact,  all 
three  key  aspects  rely  on  some  approximations  and 
include  significant  uncertainties.  Furthermore,  because 
of  its  relatively  high  frequency,  the  majority  of  the  wave 
momentum  for  an  upward-propagating  gravity  wave  is 
deposited  after  its  breakdown  (i.e.,  in  the  region  above 
its  breaking  level  Zb)- 

Alexander  and  Dunkerton  (1999)  proposed  a  simpli¬ 
fied  version  of  Lindzen’s  parameterization  scheme  by 
neglecting  the  complex  details  of  wave  dissipation  either 
below  or  above  the  breaking  level  Zb •  They  assumed  all 
momentum  of  a  wave  component  to  be  deposited  at  the 
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breaking  level  Zb-  By  doing  so,  they  mathematically  have 
directly  mapped  the  source  momentum  spectrum  in 
horizontal  wavenumber  and  phase  speed  at  the  lower 
boundary  to  a  momentum  deposition  and  eddy  diffusion 
coefficient  in  altitude.  Such  a  simplified  parameteriza¬ 
tion  not  only  increases  the  efficiency  of  computation; 
more  importantly,  it  also  combines  all  the  uncertainties 
in  the  above  three  aspects  [(i)-(iii)]  into  one  of  speci¬ 
fying  the  source  spectrum  of  the  waves.  In  principle,  this 
makes  it  much  simpler  to  validate  and  improve  the  pa¬ 
rameterization  scheme  as  more  observations  of  gravity 
wave  variance  and  spectra  become  available.  Further¬ 
more,  we  will  show  below  that  such  a  simplification  also 
makes  the  evaluation  of  the  newly  included  sensible 
heat  flux  Eq.  (lie)  straightforward.  The  3D  parame¬ 
terization  scheme  that  includes  the  breaking  trinity  in¬ 
troduced  below  essentially  combines  the  formulations  of 
Holton  and  Zhu  (1984),  Zhu  (1987),  and  Alexander  and 
Dunkerton  (1999)  to  give  a  set  of  output  profiles  of  ( Su , 
Sv,  ST)  and  ( Kzz-m ,  Kzz-T)  as  functions  of  altitude  for 
a  set  of  given  input  wind  and  temperature  profiles  of 
(u,  v,  T )  at  each  model  grid. 

From  Zhu  (1987),  the  wave  energy  density  E  for  a  3D 
hydrostatic  gravity  wave  packet  is  given  by 


E  = 


Poo^2<^2 

2  (A2-/2)2’ 


(13) 


where  p00  is  the  air  density  at  the  lower  boundary  z0  and 
O  is  the  geopotential.  Assuming  / 2  <C  co2  and  neglecting 
the  inertial  effect  in  Eq.  (13),  we  have  the  wave  action 
density  for  a  3D  gravity  wave  in  a  3D  background  flow: 


^  =  ^  =  Poo^ _ 

2K(c  —  U  cosP)3  ’ 


(14) 


where  we  have  used  Eq.  (12)  to  replace  the  wave  fre¬ 
quency  co  with  the  phase  speed  c.  We  note  that  the 
pseudomomentum  flux  FP  for  a  2D  model  defined  by 
Alexander  and  Dunkerton  (1999)  is  equivalent  to 
^/f2  +  F 2  or  KcgzA  in  the  present  notation  for  a  3D 
model  of  gravity  wave-mean  flow  interaction.  Adopting 
the  same  approximation  of  / 2  <C  co2  in  the  dispersion 
relation  (9)  and  further  assuming  A2  >  0.25 IH2  for 
breaking  waves  in  Eqs.  (9)  and  (10),  we  have 


F 


p 


=  Kc  A  = 

gz 


Poo™2 

2N(c  —  U  cosO) 


(15) 


Equation  (15)  establishes  a  close  connection  between 
the  momentum  flux  components  in  Eq.  (11)  and  the 
wave  parameters  commonly  observable  and  specified  in 
most  parameterization  schemes. 


Now,  turning  to  Holton  and  Zhu  (1984),  the  pseudo¬ 
momentum  flux  of  a  3D  wave  component  in  a  2D  flow 
[u(z),  0]  with  buoyancy  frequency  N(z)  is  given  by 


Fu  =  Pou'w'  = 


P00AOq|c  —  wcosP0| 

2 N \c  —  uQ  cos 0Q | (c  —  u  cos 0Q) 


X  exp  -2 


A  dz  , 


(16) 


where  O0  is  the  geopotential  at  the  lower  boundary  zo,  A  t 
is  the  imaginary  part  of  the  vertical  wavenumber,  and  the 
exponential  factor  in  Eq.  (16)  represents  the  damping 
effect  on  the  wave  action.  Note  that  both  A/  and  cot  rep¬ 
resent  the  dissipation  effect  of  a  wave  packet;  their  re¬ 
lationship  is  given  by  Xfi  =  — <w.A,  which  can  be  derived 
by  the  dispersion  relation  (Holton  and  Zhu  1984).  Also, 
the  wave  pseudomomentum  flux  is  inversely  proportional 
to  the  buoyancy  frequency  because  waves  are  more  easily 
amplified  in  a  less  statically  stable  atmosphere.  At  the 
lower  boundary  zo,  where  u  —  u0,  we  have 


Fp0  =  (K/k)Fu0 


Poo^o 

2 N(c  —  u0  cos 0O)  * 


(17) 


Comparing  Eq.  (17)  with  Eq.  (15)  shows  a  great  simi¬ 
larity  between  the  two  and  also  shows  how  a  formulation 
for  a  2D  flow  [z7(z),0]  can  be  appropriately  extended  to 
a  3D  flow.  It  is  the  magnitude  of  the  background  wind  U 
and  the  wavenumber  vector  direction  with  respect  to 
the  background  wind  (6  =  60  —  8)  that  matters  for¬ 
mally.  Specifically,  if  we  follow  a  procedure  similar  to 
Holton  and  Zhu  (1984)  for  a  3D  flow  [u(z),v(z)],  we  can 
derive  the  condition  for  the  determination  of  the  breaking 
level  Zb- 


F 


PO 


Po(Zft)* 

2 N(zb) 


c  —  U b  cos 6 


,3 1 c  -  UQ  cos#| 

1  (c  —  U0  cos 6) 


(18) 


where  FP0  is  given  by  Eq.  (15)  with  N  =  N0  and  U  =  £/0; 
N(zb)  and  po(zb)  are  the  background  buoyancy  fre¬ 
quency  and  air  density  at  the  breaking  level,  respec¬ 
tively.  We  have  also  assumed  a  slowly  varying  buoyancy 
frequency  in  this  derivation. 

At  this  stage,  we  follow  Alexander  and  Dunkerton 
(1999)  by  neglecting  any  damping  to  gravity  waves  be¬ 
low  the  breaking  level  [i.e.,  setting  A*  =  0  in  Eq.  (18)]. 
Then,  the  condition  for  wave  breaking  is  simplified  to 
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1/3 

exD  (Zb  ~ 

N0(c  -  UQ  cos 6) 

6XPl  3 H  ) 

(19) 


where  cb  is  the  breaking  level  intrinsic  phase  speed. 
Equation  (19)  is  a  direct  extension  of  Eq.  (6)  in  Alexander 
and  Dunkerton  (1999)  to  a  3D  wave  packet  in  a  3D  at¬ 
mosphere.  It  also  defines  the  breaking  level  intrinsic  fre¬ 
quency  a)b  =  Kcb. 

To  examine  how  the  breaking  condition  is  satisfied, 
we  rewrite  the  above  equation  at  altitude  z: 

(c  -  U cosfl)(c  -  U0cos6)m[N0/N(z)]V3 

=  <exp(^r)-  (20) 


co2  =  (co  —  ku  —  lv)2  = 


f2(  A2  +  0.25/77 2)  +  N2K2 
K2  +  A2  +  0.25/772 


(22) 


The  hydrostatic  dispersion  relation  Eq.  (9)  can  be  re¬ 
covered  from  Eq.  (22)  by  an  approximation  K2  (A2  + 
0.25/772)  for  large-scale  waves,  corresponding  to  a  hori¬ 
zontal  wavelength  much  greater  than  its  characteristic 
vertical  scale.  Equation  (22)  explicitly  shows  that  the 
squared  intrinsic  frequency  co2  is  the  weighted  average 
of  two  squared  cutoff  frequencies  / 2  and  N2.  This  leads 
to  the  following  simplified  condition  for  the  permitted 
phase  speed  for  the  default  setting  of  the  parameters  in 
the  parameterization  scheme: 


K2 


<(c-U  cos 0)2  < 


N2 

K2  +  0.25/T72  ’ 


(23) 


At  the  lower  boundary  z0 ,  the  left-hand  side  of  Eq.  (20) 
should  be  greater  than  the  right-hand  side: 

(c  ~  UQ  cosO  )413  >  .  (21) 

Wave  components  not  satisfying  the  above  condition  at 
the  lower  boundary  are  unstable  and  will  be  eliminated 
from  the  source  spectrum  in  the  parameterization  scheme 
(Alexander  and  Dunkerton  1999).  Note  that  the  left- 
hand  side  of  Eq.  (20)  represents  the  effect  of  the  critical 
level  on  wave  amplitude,  whereas  the  right-hand  side  of 
Eq.  (20)  is  the  amplification  of  the  wave  amplitude  due 
to  the  atmospheric  density  effect.  Also,  the  background 
stratification  modifies  the  wave  amplitude  so  that  the 
pseudomomentum  flux  is  inversely  proportional  to  the 
buoyancy  frequency  as  shown  in  Eq.  (16).  As  z  increases, 
the  breaking  condition  (19)  for  an  upward-propagating 
gravity  wave  in  the  absence  of  damping  (A*  =  cot  =  0)  can 
be  satisfied  through  either  the  atmospheric  density  effect 
that  increases  the  right-hand  side  of  Eq.  (20)  exponen¬ 
tially  or  the  effect  of  a  critical  level  where  the  left-hand 
side  of  Eq.  (20)  greatly  decreases  as  UcosO  — ►  c.  Equation 
(20)  also  implies  that  the  wave  action  density  is  strictly 
conserved  below  the  breaking  level  Zb-  Physically,  this 
means  that  a  wave  breaks  when  it  approaches  a  critical 
level  of  increasing  lapse  rate  because  of  the  increase  of 
the  vertical  wavenumber  and/or  when  its  amplitude 
increases  because  of  the  decrease  of  the  background  air 
density.  Furthermore,  and  most  importantly,  according 
to  Eqs.  (7)  and  (11)  the  pseudomomentum  and  sensible 
heat  fluxes  are  constant  before  the  wave  reaches  the 
breaking  level:  FP(z)  =  EP0  and  F^z(z)  =  0  for  z  <  Zb 
when  A*  =  co/  =  0. 

A  complete  dispersion  relation  for  a  nonhydrostatic 
inertia-gravity  wave  is  given  by  (Marks  and  Eckermann 
1995) 


where  the  Coriolis  parameter  /  is  a  prescribed  cutoff 
frequency  for  trapped  waves  at  the  lower  boundary  and 
N  =  N(z )  is  the  cutoff  frequency  for  reflected  waves  at 
the  entire  altitude.  Waves  that  do  not  satisfy  Eq.  (23) 
before  encountering  breaking  levels  will  also  be  elimi¬ 
nated  from  the  spectrum.  Condition  (23)  is  again  a  3D 
extension  of  a  similar  condition  of  total  internal  re¬ 
flection  when  co2  exceeds  N2  given  by  Alexander  and 
Dunkerton  (1999)  and  also  the  possible  reflection  by 
Jones’  critical  levels  of  co2  =  f 2  (Yamanaka  and  Tanaka 
1984)  at  the  lower  boundary.  Computationally,  we  fur¬ 
ther  set|/|>2XlO-5  s_1  so  that  the  left-side  inequality 
of  Eq.  (23)  excludes  unusually  large  values  of  Ffo  in 
specifying  the  source  spectra  based  on  Eq.  (15)  for  a 
given  O0  field  at  the  lower  boundary. 

For  waves  that  reach  breaking  levels,  we  mainly  fol¬ 
low  Alexander  and  Dunkerton  (1999)  and  assume  that 
all  the  momentum  is  deposited  at  the  breaking  levels, 
producing  the  breaking  trinity  of  momentum  drag,  eddy 
diffusion,  and  wave  heating  locally.  Mathematically,  the 
procedure  is  equivalent  to  mapping  the  source  mo¬ 
mentum  spectrum  in  horizontal  wavenumber  and  phase 
speed  to  a  momentum  deposition  in  altitude.  Such  a  map¬ 
ping  procedure  also  makes  the  evaluation  of  Eq.  (lie) 
straightforward.  To  illustrate  this  point  explicitly,  we  re¬ 
write  Eq.  (lie)  as 

f*.*  -  ■  -nD^Kc<«A)- 

\  ZZ  gzj 

(24) 

where  Dzz  is  the  wave  eddy  diffusion  coefficient  at  the 
breaking  level  Zb  and  is  proportional  to  the  momen¬ 
tum  deposition  rate  for  the  breaking  wave  component 
(Holton  1982;  Alexander  and  Dunkerton  1999): 
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Dzz  =  -  U  cosd)/N2.  (25) 

Below  the  breaking  level,  F ^  z  vanishes  because  (ot  =  0 
for  a  wave  packet  without  dissipation.  At  z  =  Zb,  <*>i  = 
A 2Dzz  (Lindzen  1981)  so  we  have  a  nonzero  heat  flux 
once  the  wave  starts  breaking.  Note  that  for  a  given 
horizontal  wavenumber  K ,  the  intrinsic  phase  speed  cb 
and  frequency  cob  at  the  breaking  level  are  known  and 
determined  by  Eq.  (19).  As  a  result,  the  corresponding  A 
and  cgz  at  Zb  can  also  be  easily  calculated  by  the  dis¬ 
persion  relations  [Eqs.  (9)  and  (10)]: 


terms  in  Eqs.  (lla)-(llc)  can  be  simultaneously  and 
self-consistently  calculated  once  the  pseudomomentum 
flux  KcgzA  is  specified  at  the  lower  boundary  for  a  set  of 
wave  parameters.  In  the  current  notation,  the  induced 
breaking  trinity  of  momentum  drag,  eddy  diffusion  co¬ 
efficient,  and  wave  heating  rate  at  altitude  Zk  by  break¬ 
ing  wave  components  (/  =  1, . . . ,  J)  can  be  calculated  in 
the  following  formulations: 

Xh  =  pJ7jAz^(CO&fti]Fm)i  for  z^<zh{ci)  - 

(27) 


2  N2 

A  **= 


K2N2 


A3<w, 


(26a, b)  Y  _  e 

h  Po(zh)Az 


X(sineoFpo)y  for  zk_x  <zb{cj)  <  zk. 


Therefore,  the  coefficient  defined  in  brackets  in  Eq. 
(24)  can  be  expressed  in  known  wave  parameters: 
II  =  N3I(c\K).  In  Eq.  (26),  we  have  again  assumed  A  to 
be  much  greater  than  1/(2 H)  at  Zb-  This  is  a  good  ap¬ 
proximation  because  waves  generally  break  near  critical 
levels  where  A2  >  (0.25/H2).  At  this  stage,  all  the  flux 


Dh  = 


Po(h)N  ( zh)Az 


X[c  -  U(zh)cosd]j(Fpo)j 


for 


Zk-l  <  zb(cj )  -  zv 


(28) 


(29) 


Qh 


Y 

P2o(zh)N2(zh)Az2  i 

y 

.Po(zh)N2(zh)Az2  i 


[c-U(zh)cos0]j(UF2o)j, 

\c-U(zh)cos6]i(HF2n)r 


zk-i<z  ~  zk 


zk<z~  zk+ 1 


(30) 


where  Zh  =  Zk  ~  Az/ 2  and  Zh-i  =  Zk  ~  Az  represent  the 
half-  and  full-grid  step  below  Zk ,  respectively.  Outside  the 
breaking  regions  as  specified  in  the  above  expressions, 
the  values  are  set  to  zero  for  those  breaking  waves.  The 
force  terms  on  the  right-hand  sides  of  Eqs.  (l)-(3)  on  the 
model  grids  are  the  averages  of  the  above  expressions: 


(31) 

+  U  i). 

(32) 

Kzz-m(zk)  =  2^h  +  Dh+i)’ 

(33) 

Kzz_T(zk)  =  Kzz_JZk)/Pr, 

(34) 

S Azk )  =  2^h  +  ®h+l)- 

(35) 

In  the  above  equations,  s  is  an  intermittency  factor 
introduced  by  Alexander  and  Dunkerton  (1999)  that 
represents  the  ratio  of  the  observed  momentum  flux  to 


the  modeled  one  at  the  lower  boundary.  It  is  a  specified 
parameter  in  the  parameterization  scheme  that  can  vary 
with  time,  space,  and  wave  parameters.  The  eddy  dif¬ 
fusion  coefficient  derived  from  the  Lindzen-type  pa¬ 
rameterization  refers  to  the  dissipation  and  momentum 
mixing  to  the  wave  field.  To  apply  it  to  the  background 
mean  state,  we  introduce  an  additional  parameter  sm  to 
characterize  its  efficiency  for  mixing  the  momentum  of 
the  background  fields.  Holton  and  Zhu  (1984)  showed 
that  there  was  a  partial  cancellation  between  the  direct 
wave  drag  and  the  drag  induced  by  the  eddy  diffusion, 
which  smoothed  the  total  drag  on  the  right-hand  sides  of 
Eqs.  (l)-(3).  The  current  parameterization  scheme  in¬ 
cludes  hundreds  to  thousands  of  wave  components  in  the 
source  spectrum,  which  yields  smooth  drag  profiles  for 
typical  wind  profiles.  Therefore,  we  set  a  smaller  value  of 
sm  =  0.3  in  the  current  scheme.  The  eddy  Prandtl  number 
Pr  is  defined  as  the  ratio  of  the  eddy  momentum  diffu- 
sivity  to  the  eddy  heat  and  tracer  diffusivity.  It  also  reflects 
the  fact  that  the  wave  eddy  diffusion  coefficient  Dzz  defined 
in  Eq.  (24)  is  different  than  the  eddy  diffusion  coefficient 
for  the  background  state  Kzz^m  defined  in  Eq.  (33).  Both 
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the  mesospheric  tracer  measurements  and  careful  model 
analysis  of  the  breaking  process  show  that  Pr  is  much 
greater  than  1  (e.g.,  Strobel  et  al.  1987;  Strobel  1989;  Chao 
and  Schoeberl  1984).  This  is  consistent  with  Eq.  (8)  of  a 
vanishing  sensible  heat  flux  to  a  first-order  approximation. 
We  set  Pr  =  5  in  the  current  scheme.  Again,  it  is  noted 
from  Eqs.  (27)  and  (28)  that  the  partition  of  the  deposited 
momentum  in  (x,  y)  direction  is  based  neither  on  (cos5, 
sinS)  nor  on  (cos#,  sin 0)  but  rather  on  (cos0o,  sinflo)* 

It  should  be  pointed  out  that,  strictly  speaking,  the  first 
and  second  terms  on  the  right-hand  sides  of  Eqs.  (l)-(3) 
represent  two  different  ways  of  parameterizing  the  ef¬ 
fects  of  subgrid-scale  motions:  nonlocalized  wave  forcing 
due  to  wave  propagation  and  breaking  versus  localized 
mixing-length  theory  of  turbulence.  For  a  single  gravity 
wave  component,  the  wave-breaking  processes  occurring 
on  scales  much  smaller  than  its  wavelength  dissipate  the 
wave  action  and  induce  the  explicit  momentum  and  heat 
flux  divergences  shown  in  Eqs.  (4)  and  (11).  When 
there  is  a  clear  scale  separation  among  the  background 
grid-resolved  motions,  the  smaller  subgrid  wave  compo¬ 
nent,  and  a  much  smaller-scale  mixing  by  wave  breaking, 
the  effect  of  the  subgrid-scale  eddy  terms  can  be  entirely 
included  by  the  first  terms  (Su,  Sv,  ST )  in  Eqs.  (l)-(3). 
Including  the  second  eddy  diffusion  terms  on  the  right- 
hand  sides  of  Eqs.  (l)-(3)  implies  that  the  vertical  mixing 
by  the  much  smaller-scale  motions  that  directly  dissipates 
the  wave  component  also  diffusively  mixes  the  grid- 
resolved  motions.  The  diffusion  terms  are  expected  to 
become  more  important  as  more  gravity  waves  compo¬ 
nents  with  different  wavelengths  are  included  in  the 
parameterization  scheme  because  this  will  make  the 
separation  among  different  scales  less  obvious.  This 
analysis  also  suggests  that  the  additional  parameter  sm 
introduced  in  the  current  parameterization  scheme  is 
empirical  and  it  may  vary  with  the  setting  of  the  wave 
source  spectrum. 

From  Eqs.  (15),  (19),  (21),  and  (26),  we  see  that,  at  Zb , 
c  and  thus  co ,  EP0,  and  II  are  all  positive.  Therefore,  the 
signs  of  the  drag  components  are  solely  determined  by 
the  wavenumber  vector  direction  (cos0o,  sin0o),  whereas 
the  heating  rate  and  eddy  diffusion  coefficients  are  al¬ 
ways  positive.  Previous  studies  (e.g.,  Walterscheid  1981; 
Liu  2000;  Talaat  et  al.  2001)  suggested  that  the  gravity 
waves  could  induce  dynamical  cooling  in  the  wave- 
dissipating  region.  A  positive  heating  at  Zb  in  the  current 
parameterization  scheme  is  mainly  due  to  the  assumption 
that  each  spectral  wave  component  is  completely  dissi¬ 
pated  at  the  breaking  level  where  the  wave  heating  occurs. 
On  the  other  hand,  if  one  assumes  a  finite  and  a  broad 
region  of  wave  dissipation  between  Zb  and  a  higher  level 
such  as  a  critical  level  zc  where  the  wave  amplitude 
completely  vanishes,  then  cooling  could  occur  in  most 


of  the  region  between  Zb  and  zc-  To  accommodate  such  a 
cooling  effect,  we  add  a  cooling  term  of  the  same  mag¬ 
nitude  in  Eq.  (30)  immediately  above  Zb,  which  leads  to 
a  net  effect  of  downward  transport  of  heat  near  the  wave 
breaking  level.  Physically,  this  means  that  there  is  no  net 
heat  being  induced  within  the  atmosphere  when  its 
vertical  heat  flux  vanishes  in  both  the  lower  and  upper 
boundaries.  Mathematically,  the  heat  flux  term  of 
an  upward-propagating  gravity  wave  packet  shown  in 
Eq.  (lie)  experiences  two  major  stages:  (i)  it  changes 
from  0  to  a  finite  value  as  the  wave  breaks  because  of 
a  finite  jump  in  (ot  at  Zb  and  (ii)  it  changes  back  to  0  again 
as  the  wave  approaches  a  critical  level  and  completely 
dissipates  so  that  coA  —>  0.  Previous  studies  mainly  fo¬ 
cused  on  (ii)  and  overlooked  (i).  The  net  effect  of  Eq.  (30), 
where  the  breaking  waves  heat  the  background  atmosphere 
at  Zb  and  cool  it  one  level  above  Zb,  is  the  downward 
transport  of  heat.  We  point  out  that  the  physical  mech¬ 
anism  of  this  downward  heat  transport  caused  by  a  dis¬ 
continuous  wave  breaking  is  slightly  different  from  the 
same  conclusion  in  previous  studies  (e.g.,  Walterscheid 
1981),  where  the  downward  heat  flux  and  the  associated 
cooling  only  correspond  to  the  above  stage  (ii). 

Alexander  and  Dunkerton  (1999)  emphasized  the  dy¬ 
namical  importance  of  wave  dissipation  at  Zb  rather  than 
at  zc-  The  above  analysis  of  the  heat  flux  and  transport 
shows  the  energetic  importance  of  wave  dissipation  at  Zb- 
Wave  dissipation  at  Zb  produces  a  net  downward  trans¬ 
port  of  thermal  energy  to  the  mean  state.  Realization  of 
mi  =  0  and  coA  =  0  at  the  lower  and  upper  boundaries 
respectively  allows  the  finite  downward  heat  transport  to 
be  evaluated  self -consistently.  This  also  makes  Eq.  (3) 
energetically  consistent.  In  summary,  the  current  pa¬ 
rameterization  of  the  breaking  trinity  extends  the  pre¬ 
vious  Lindzen  types  of  parameterizations  such  as  those  by 
Holton  and  Zhu  (1984)  and  Alexander  and  Dunkerton 
(1999)  in  two  aspects:  extending  it  to  a  3D  background 
flow  and  including  a  wave  breaking-induced  heating 
term.  For  a  given  set  of  input  wind  and  temperature 
profiles  [u(z),  v(z),  T(z)\  as  a  basic  state  in  a  model  grid, 
the  parameterization  scheme  of  the  breaking  trinity 
outputs  the  vertical  profiles  of  (Su,  Sv ,  ST,  Kzz-m ,  Kzz-T ) 
that  allow  one  to  calculate  all  the  force  terms  by  the 
subgrid-scale  motions  such  as  those  on  the  right-hand 
sides  of  Eqs.  (l)-(3). 

c.  Specification  of  the  source  spectrum 

A  discrete  source  spectrum  needs  to  be  specified  for  the 
geopotential  variance  Oo(0o,  c,  K)  as  a  function  of  three 
wave  parameters:  0O,  c  and  K.  We  have  already  indicated 
the  reason  for  choosing  these  three  independent  wave 
parameters  for  prescribing  the  wave  spectrum:  0O,  c,  and 
K  do  not  vary  with  altitude  for  an  upward-propagating 
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wave  component.  In  some  parameterization  schemes, 
such  as  the  one  by  Warner  and  McIntyre  (1996),  the 
distribution  and  variation  of  the  energy  spectrum  were 
prescribed  as  a  function  of  A  and  co,  both  of  which  vary 
with  the  altitude.  We  assume  a  functional  form  that  is 
separated  in  wavenumber  vector  direction:  Oo(0o,  c,  K)  = 
B1(9q)B2(c,  K ).  Following  Matsuno  (1982)  and  Holton 
and  Zhu  (1984),  we  assume  an  isotropic  distribution  of  Oq 
in  0o  [i.e.,  51(0Q)  =  B  =  Constant].  Depending  on  spe¬ 
cific  model  experiments  and  applications,  the  wave- 
number  vector  azimuthal  angle  60  is  divided  into  12, 16,  or 
20  equal  sectors.  In  addition  to  simplicity  (eliminating 
one  set  of  the  tuning  parameters  in  the  parameterization 
scheme),  there  are  several  physical  justifications  for  an 
isotropic  spectrum:  (i)  waves  with  small  phase  speed  with 
random  background  wind  will  climatologically  dominate 
the  topographic  gravity  waves  at  the  lower  boundary; 

(ii)  waves  induced  by  convection  (e.g.,  Alexander  et  al. 
1995)  or  by  geostrophic  adjustment  forced  by  a  local 
momentum  source  (e.g.,  Zhu  and  Holton  1987)  show  near 
isotropic  distribution  in  perturbation  geopotential;  and 
most  importantly  (iii)  the  gravity  wave  source  spectrum 
added  to  the  model  will  not  change  its  dynamical  prop¬ 
erty  of  conserving  the  total  angular  momentum.  Note 
from  Eq.  (15)  that  there  exists  a  difference  in  definition  of 
the  isotropic  spectrum  between  the  geopotential  variance 
Oq  and  momentum  flux  FP.  However,  the  first  two  justi¬ 
fications  described  above  (i.e.,  small  c  with  random  U  and 
large  c  with  small  U)  also  lead  to  an  isotropic  distri¬ 
bution  in  FP . 

We  assume  a  wide  range  of  horizontal  wavenumber 
with  12  equally  spaced  sectors  in  logarithmic  scale  be¬ 
tween  27r/(800  km)  and  27t/(6.25  km).  This  logarithmic 
range  is  the  same  as  the  one  used  by  Alexander  (1998) 
in  a  modeling  study  of  linear  wave  ray  tracing  that  re¬ 
produces  several  climatological  patterns  observed  for 
stratospheric  gravity  wave  variance.  Such  a  setting  in¬ 
corporates  several  well-known  observational  and  mod¬ 
eling  facts  about  gravity  waves:  (i)  the  characteristic 
wavelengths  of  the  terrain-generated  gravity  waves  are 
on  the  order  of  a  few  to  tens  of  kilometers  (e.g.,  Nappo 
2002);  (ii)  the  dominant  wavelengths  of  the  convectively 
generated  gravity  vary  with  period  and  range  from 
tens  to  hundreds  of  kilometers  (e.g.,  Alexander  1996); 

(iii)  the  typical  horizontal  wavelength  generated  by  geo¬ 
strophic  adjustment  under  a  localized  momentum  forcing 
is  about  300  km  (Zhu  and  Holton  1987);  and  (iv)  in¬ 
creasing  the  horizontal  wavelength  reduces  Zb  (Alexander 
and  Dunkerton  1999),  which  allows  for  a  deeper  and  more 
uniform  distribution  of  the  parameterized  drag  and  eddy 
diffusion.  We  also  note  that  the  typical  horizontal  res¬ 
olution  of  current  GCMs  is  on  the  order  of  hundreds 
of  kilometers.  The  specification  of  the  phase  speed 


spectrum  B2(c ,  K)  follows  Matsuno  (1982),  Alexander 
and  Dunkerton  (1999),  and  most  modeling  studies  that 
use  a  Gaussian  function  to  construct  the  spectrum: 


BJc,  K)  =  X  W.  exp 


(36) 


where  the  central  phase  speed  c0j,  half -width  cwj,  and  the 
corresponding  weight  Wj  for  different  horizontal  wave- 
number  Kj  (j  =  1,  . . .  12)  are  specified  according  to  var¬ 
ious  observational  and  modeling  constraints.  Because  the 
direction  of  the  wave  propagation  has  been  represented 
by  the  wavenumber  vector  direction  0O,  both  c  and  K  are 
positive.  In  the  default  setting,  the  phase  speed  is  divided 
equally  into  16  sectors  between  5  and  70  m  s-1,  with  c0,7 
and  cwj  varying  linearly  within  [0,  50]  and  [10,  20],  re¬ 
spectively.  Computationally,  the  default  setting  of  the 
three  wave  parameters  consists  of  20  equal  sectors  in  0O, 
16  intervals  in  c,  and  12  equally  spaced  sectors  in  loga¬ 
rithmic  scale  in  K.  Therefore,  the  source  spectrum  con¬ 
tains  3840  waves  (=20  X  16  X  12).  Most  of  these  waves 
will  break  at  different  altitude  levels  to  produce  smooth 
profiles  for  the  source  terms  in  Eqs.  (l)-(3). 

The  entire  spectrum  is  normalized  by  the  measure¬ 
ments  based  on  the  total  variance  of  the  gravity  wave 
geopotential  Oq  at  the  lower  boundary,  which  varies 
with  space  and  time  and  is  constrained  by  the  observa¬ 
tional  climatology.  To  give  an  estimate  of  its  magnitude, 
we  note  that  the  variance  of  the  geopotential  height  for 
planetary  waves  at  250-hPa  in  mid-  and  high-latitude 
regions  is  about  100-200  m  (e.g.,  James  1994).  If  we 
assume  that  the  variance  of  the  geopotential  height  for 
gravity  waves  at  the  lower  boundary  is  two  orders  of 
magnitude  smaller  than  that  of  the  planetary  waves  at 
the  troposphere,  then  Oq  is  on  the  order  of  (10)2 
(m2  s-2)2  or  (20)2  (m2  s-2)2.  In  Alexander  and  Dunkerton 
(1999),  it  is  suggested  that  over  topography,  the  mag¬ 
nitude  of  the  observed  pseudomomentum  flux  p0u'w'  is 
0.03  -  0.5  Pa.  From  Eq.  (15),  the  variance  of  the  cor¬ 
responding  geopotential  is 


O2 


2N\c  —  UQ  cos0| 
PooK 


P0 


'  500  (m2  s"2)2  (37) 


for  a  typical  setting  of  N  =  0.02  s-1,  p0 0  =  1.29  kg  m-3 
at  the  surface,  | c  -  U0  cos0|  =  20  m  s-1,  and  K  =  2ttI 
(100  km),  which  is  consistent  with  the  above  estimate  scaled 
according  to  the  planetary  wave  activity  in  the  troposphere. 


3.  Numerical  results  of  the  parameterization 
scheme 

A  set  of  output  fields  near  the  Northern  Hemisphere 
solstice  (20050105)  from  version  5  of  the  Goddard  Earth 
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Fig.  1.  Example  of  the  winds  from  a  high-altitude  version  of  the  GEOS-5  model  (20050105)  and  the  parameterized  breaking  trinity: 
wave  drag,  eddy  diffusion  coefficient,  and  wave  heating  by  gravity  wave  breaking  from  the  GEOS-5  model  runs  near  the  NH  winter 
solstice.  All  fields  have  been  zonally  averaged  over  144  longitudinal  grids:  (top  left)  zonal  and  (top  right)  meridional  wind;  (middle  left) 
zonal  and  (middle  right)  meridional  drag;  and  (bottom  left)  eddy  diffusion  coefficient  and  (bottom  right)  heating  rate. 


Observing  System  Data  Assimilation  System  (GEOS-5 
DAS)  is  used  to  test  the  parameterization  scheme.  The 
GEOS-5  atmospheric  GCM  is  a  weather-climate-capable 
model  consisting  of  a  finite-volume  dynamical  core  and 
a  physical  package  parameterizing  four  major  groups  of 
physical  processes  (Rienecker  et  al.  2008).  The  standard 
setting  of  the  atmospheric  model  has  72  vertical  layers 
extending  from  the  surface  to  —70  km.  The  major  physical 


processes  contained  in  the  model  are  (i)  moist  processes 
including  cloud  microphysics,  (ii)  shortwave  and  infrared 
radiation,  (iii)  drag  and  eddy  diffusion  parameterized  by 
a  2D  Lindzen-type  of  scheme  for  gravity  wave  breaking 
(McFarlane  1987;  Garcia  1991),  and  (iv)  surface  pro¬ 
cesses  in  the  atmospheric  boundary  layer.  A  high-altitude 
version  of  the  GEOS-5  atmospheric  model  has  82  ver¬ 
tical  levels  extending  from  the  surface  to  —100  km  with 
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Fig.  2.  Example  of  the  parameterized  breaking  trinity:  wave  drag,  eddy  diffusion  coefficient,  and  wave  heating  by  gravity  wave  breaking 
derived  from  the  zonally  averaged  wind  and  temperature  fields  near  the  NH  winter  solstice:  (top  left)  zonal  and  (top  right)  meridional 
drag;  and  (bottom  left)  eddy  diffusion  coefficient  and  (bottom  right)  heating  rate. 


a  horizontal  resolution  of  2.5°  longitude  by  2°  latitude.  In 
Fig.  1,  we  show  the  zonal  mean  fields  of  the  assimilated 
model  winds  and  the  parameterized  breaking  trinity 
calculated  by  the  current  scheme  based  on  the  output 
model  winds:  wave  drag,  eddy  diffusion  coefficient,  and 
wave  heating  by  gravity  wave  breaking.  We  do  not  show 
the  temperature  field  in  the  figure  because  it  only  plays 
a  minor  role  in  the  parameterization  scheme  by  modify¬ 
ing  the  background  static  stability.  Note  that  the  magni¬ 
tudes  of  the  zonal  mean  zonal  drag  Su  and  eddy  diffusion 
coefficient  Kzz_m  around  the  mesopause  are  comparable 
to  those  in  previous  studies  based  on  traditional  Lindzen- 
type  gravity  wave  parameterizations,  which  were  re¬ 
quired  to  produce  the  horizontal  anomalous  temperature 
gradient  (or  the  vertical  zonal  wind  reversals)  near  the 
solstice  mesopause  (e.g.,  Holton  1983).  In  addition, 
the  current  3D  parameterization  scheme  also  provides 
the  forcing  components  of  the  meridional  wave  drag  Sv 
and  wave  heating  rate  ST  in  the  primitive  Eqs.  (l)-(3). 
Near  the  mesopause,  the  overall  magnitude  of  the  wave 
heating  and  cooling  due  to  wave  breaking  is  comparable  to 
other  estimates  based  on  energy  dissipation  rates  either 


from  the  model  (e.g.,  Liu  2000;  Becker  and  Schmitz  2002) 
or  from  the  measurements  (e.g.,  Liibken  1997)  and  is  also 
comparable  to  the  radiative  and  chemical  heating  rates 
in  the  same  region  (Zhu  1994;  Zhu  et  al.  2000).  In  Liu 
(2000),  the  vertical  wave  heating  rate  distribution  also 
shows  some  cancellation  in  altitude.  Note  that  the  added 
wave  heating  ST  not  only  directly  changes  the  thermal 
structure  near  the  mesopause  but  also  modifies  the  me¬ 
ridional  circulation  through  its  thermal  drive  of  the  me¬ 
ridional  gradient  of  the  heating  rate  (e.g.,  Zhu  et  al.  2001). 
Historically,  parameterizations  of  Su  and  Kzz_m  by  wave 
breaking  were  proposed  to  explain  the  observed  large- 
scale  features  of  the  horizontal  anomalous  temperature 
gradient  (or  the  vertical  zonal  wind  reversals)  near  the 
solstice  mesopause.  Since  Sv  and  ST  have  been  self- 
consistently  derived  together  with  Su  and  Kzzm ,  their 
additional  effects  on  middle  atmosphere  dynamics  and 
physics  should  also  be  an  important  subject  to  be  care¬ 
fully  examined.  It  should  be  emphasized  that  all  terms  Su , 
Sv,  Kzz_m ,  and  ST  have  been  consistently  calculated  on 
each  model  grid  simultaneously  for  the  corresponding 
individual  input  wind  and  temperature  profiles.  The  2D 
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Fig.  3.  As  in  Fig.  1,  but  on  a  longitudinal  plane  of  110°E,  where  the  meridional  wind  is  peaked. 


plots  shown  in  Fig.  1  represent  the  zonal  mean  fields  av¬ 
eraged  over  144  longitudinal  grids. 

Figure  2  shows  the  breaking  trinity  of  Su,  Sv,  Kzz_m, 
and  ST  as  a  function  of  latitude  and  altitude  calculated 
based  on  the  zonally  averaged  wind  field  as  shown  on  the 
top  row  of  Fig.  1  and  the  corresponding  averaged  tem¬ 
perature  field.  Note  that  the  wave-mean  flow  inter¬ 
actions  as  described  in  the  last  section  by  Eqs.  (19),  (23), 
and  (26)-(30)  are  highly  nonlinear.  Therefore,  the  con¬ 
sequences  of  the  breaking  trinity  calculated  based  on  a 


zonally  averaged  flow  will  not  be  the  same  as  those  shown 
in  Fig.  1,  which  were  calculated  by  zonally  averaging  the 
derived  force  terms.  Specifically,  the  overall  magnitudes 
in  Fig.  2  are  significantly  greater  than  those  in  Fig.  1. 
However,  a  careful  comparison  of  the  two  figures  shows 
that  they  have  similar  spatial  patterns.  For  example,  near 
the  solstice  mesopause,  the  zonal  momentum  drag  is 
mostly  easterly  in  the  winter  hemisphere  and  westerly  in 
the  summer  hemisphere,  which  is  dynamically  required  to 
produce  the  horizontal  anomalous  temperature  gradient 
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near  the  solstice  mesopause.  In  addition,  distribution  of 
the  wave  heating  ST  shows  well-organized  dipole  pat¬ 
terns  with  cooling  regions  sitting  above  the  heating  regions 
in  both  figures,  as  expected  from  its  parameterization 
scheme  (30).  Physically,  the  existence  of  these  differ¬ 
ences  and  commonalities  between  the  two  sets  of  mean 
fields  is  not  surprising  because  the  zonal  mean  state  in 
a  3D  numerical  model  is  driven  somewhat  differently 
than  the  one  in  a  2D  model  since  they  have  entirely 
different  subgrid  scales. 

The  wind  distribution  at  a  particular  longitude  sector 
could  significantly  deviate  from  its  zonal  average.  This  is 
similar  to  the  case  in  which  a  measured  single  wind 
profile  at  a  prescribed  location  could  be  completely  dif¬ 
ferent  than  its  climatology.  To  see  such  an  effect  and 
further  illustrate  the  3D  nature  of  the  current  parame¬ 
terization  scheme,  we  show  in  Fig.  3  the  input  wind  field 
(u,  v)  and  the  output  force  terms  ( Su ,  Sv,  Kzz_m,  ST )  on 
a  longitudinal  plane  (110°E)  where  the  meridional  wind 
is  peaked.  Comparison  of  Figs.  1  and  3  shows  that  both 
the  input  and  output  fields  have  significant  differences  in 
their  magnitudes  and  patterns.  Specifically,  the  zonal 
winds  in  the  winter  hemisphere  around  the  stratopause 
are  of  opposite  signs.  As  a  result,  the  zonal  drag  force 
near  the  winter  mesopause  also  changes  sign  because  of 
the  filtering  effect  of  the  upward-propagating  gravity 
waves.  Since  the  wave  filtering  effect  is  believed  to  be 
one  of  the  main  driving  mechanisms  for  a  variety  of 
middle  atmosphere  circulations  (e.g.,  Dunkerton  1982; 
Smith  1996;  Mayr  et  al.  1998),  the  spectral  parameteri¬ 
zation  of  wave  breaking  proposed  in  the  current  paper  is 
expected  to  be  appropriate  for  simulations  of  various 
middle  atmospheric  circulations. 

4.  Conclusions 

In  this  paper,  we  have  developed  a  3D  spectral  pa¬ 
rameterization  scheme  to  self-consistently  include  the 
“breaking  trinity”  of  upward-propagating  gravity  waves 
for  large-scale  numerical  models:  (i)  momentum  drag 
that  represents  the  nonlocalized  transport  of  momen¬ 
tum  through  wave  propagation  in  a  3D  background  flow, 

(ii)  eddy  diffusion  coefficients  that  characterize  the  lo¬ 
calized  diffusive  transport  of  momentum,  heat,  and 
tracers  due  to  3D  mixing  induced  by  wave  breaking,  and 

(iii)  wave  heating  rate  that  captures  localized  transport 
of  heat  by  perturbing  wave  structures  to  redistribute  the 
thermal  energy  within  a  finite  domain.  For  a  set  of  given 
input  wind  and  temperature  profiles  at  each  model  grid 
(u,  v ,  T ),  the  parameterization  scheme  returns  five  ver¬ 
tical  profiles  Su,  Sv,  Kvv  Kvy _T,  and  ST  for  calculating 
the  force  terms  on  the  right-hand  sides  of  the  momen¬ 
tum  and  energy  Eqs.  (l)-(3). 


The  spectral  parameterization  has  been  developed  by 
using  a  general  relationship  between  the  wave  action 
flux  and  the  wave  momentum  and  heat  fluxes  developed 
by  Zhu  (1987)  and  a  mapping  approximation  between 
the  wave  source  spectrum  and  the  vertical  distribution  of 
the  momentum  deposition  developed  by  Alexander  and 
Dunkerton  (1999).  When  the  parameterization  algo¬ 
rithm  is  applied  to  a  set  of  3D  wind  fields  output  from 
a  high-altitude  version  of  the  GEOS-5  atmospheric 
model,  the  derived  zonal  mean  drag  and  eddy  diffusion 
coefficient  near  the  solstice  mesopause  are  comparable 
to  those  derived  in  previous  work  and  are  required  to 
produce  the  horizontal  anomalous  temperature  gradient 
and  the  vertical  zonal  wind  reversals  near  the  solstice 
mesopause  (e.g.,  Holton  1983).  The  derived  wave  heating 
and  cooling  rates  near  the  mesopause  are  found  to  be  10- 
20  K  day-1,  which  is  comparable  to  the  radiative  and 
chemical  heating  rates.  The  filtering  effect  of  the  upward- 
propagating  gravity  waves  is  also  well  captured  by  the 
current  parameterization  scheme. 
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